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Backward waves and negative refraction are shown to exist in plasmonic crystals whose lattice cell 
size is a very small fraction of the vacuum wavelength (less than l/40th in an illustrative example). 
Such "quasi-homogeneity" is important, in particular, for high-resolution imaging. Real and complex 
Bloch bands are computed using the recently developed finite-difference calculus of "Flexible Local 
Approximation MEthods" (FLAME) that produces linear eigenproblems, as opposed to quadratic 
or nonlinear ones typical for other techniques. FLAME dramatically improves the accuracy by 
incorporating local analytical approximations of the solution into the numerical scheme. 



Backward waves (Poynting vector opposite to phase 
velocity [29] ) and the closely related phenomenon of neg- 
ative refraction have been extensively studied in recent 
years (e.g. [9, 12] and references there) due to the intrigu- 
ing physical effects and potential applications in imaging 
and other areas. 

Backward waves in periodic structures may in general 
exist only if the lattice cell size, as a fraction of the vac- 
uum wavelength Aq, is above certain thresholds derived 
recently in [23]. Plasmonic crystals [4, 14] are an in- 
teresting exception because in the vicinity of a plasmon 
resonance the constraints on the cell size are relaxed or 
removed, as explained below, and thus it may be possi- 
ble to reduce the lattice cell size and approach an ideal 
homogeneous negative-index medium. 

Bloch modes play a central role in the analysis of elec- 
tromagnetic waves in periodic structures. While ana- 
lytical expressions for these modes are available only in 
special one-dimensional cases [11, 22, 26], there exist a 
variety of computational techniques: Fourier transforms 
(plane wave expansion, PWE) [7, 11, 22], scattering ma- 
trices and lattice summation [2], finite difference [25, 27] 
and finite element [1, 4, 5, 8, 17] analysis, semi-analytical 
methods [28], and more. 

Bloch wave problems have three (in 2D) or four (in 
3D) scalar eigenparameters: frequency uj and the Carte- 
sian components of the Bloch vector K. Solving for all 
these parameters, and the respective eigenmodes, simul- 
taneously is impractical. The usual approach is to look 
for the values of lo for any given K. The differential op- 
erator of the problem, and hence the respective matrices 
in the numerical computation, contain the permittivity 
e and therefore for dispersive media (e = e('j^)) depend 
on the frequency in a complicated way. Consequently, 
the resulting eigenvalue problems with respect to lo are 
nonlinear. Several solution methods have been proposed 
[16, 19] but are not simple, and convergence is not guar- 
anteed. 

A more elegant approach, where frequency is treated 
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as a given parameter and components of the Bloch vec- 
tor as unknown eigenvalues, has been explored relatively 
recently in PWE [13] and in FEM [4, 17]. Quadratic 
eigenproblems with respect to the Bloch number usually 
arise and can be converted to linear ones by introduc- 
ing auxiliary unknowns either on the continuous level 
(e.g. solving for both fields E, H instead of just one) 
or, alternatively, on the linear algebra level [18]. This 
conversion doubles the number of unknowns; the compu- 
tational cost, typically proportional to the cube of the 
system size, increases about eightfold. 

In the recently developed generalized finite-difference 
(FD) calculus of Flexible Local Approximation MEthods 
(FLAME [21, 22, 24]) high accuracy is achieved by re- 
placing the Taylor expansions of standard FD analysis 
with much better approximating functions, e.g. plane 
waves or cylindrical harmonics. In FLAME, a; is a nat- 
ural "independent variable" because the approximating 
functions are derived for a fixed value of w. FLAME has 
two clear advantages in the computation of (real or com- 
plex) Bloch modes: (i) it dramatically improves the ac- 
curacy by incorporating local analytical approximations 
of the solution into the numerical scheme (see examples 
below); (ii) it produces linear eigenproblems. 

Let us consider band structure calculation in a pho- 
tonic crystal formed by an infinite lattice of rectangu- 
lar cells Lx X Ly in the xy-plane. In a very common 
case, each cell contains a dielectric cylindrical rod with 
a radius Trod and the relative dielectric permittivity Crod- 
The medium outside the rod has permittivity tout- At 
optical frequencies, all media are intrinsically nonmag- 
netic. The governing wave equation for the TE mode 
(one-component magnetic field phasor H = H^) is 

V-e-^VH + Lo^fioH = (1) 

where the relative dielectric permittivity e — e{x,y) is 
periodic over the lattice. The exp(—iujt) convention is 
used for complex phasors. The H field is sought as a 
Bloch-Floquet wave [11] with a (yet undetermined) Bloch 
vector Kb- 

H{r) = i/pER(r) exp(iKs • r) r = {x,y) (2) 

where -ffpER is periodic over the lattice. In the 
space of Bloch vectors K^, the first Brillouin zone is 
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FIG. 1: A fragment of the band diagram of the Davanco et 
al. [4] plasmonic crystal. Cu — ua/c, K = KBa/T^- 



[—'k/Lx,t^/Lx\ X [— 7r/Lj^, tt/Lj,]. For notational simplic- 
ity and without real loss of generality, let L^ ^ Ly = a. 

There are two general options: solving for the periodic 
factor iJpER(a;,y) or, alternatively, for the full iZ-field of 
(1). In the first case, standard periodic boundary condi- 
tions apply, but the differential operator is more compli- 
cated than in the second case. The boundary conditions 
for the full iJ-field are "scaled-periodic" due to the Bloch 
exponential exp(iK • r): 

H{a/2,y)^exp{iK,a)H{-a/2,y); \y\<a/2 (3) 
with similar conditions at the boundaries y — ±a/2. 

Accurate local analytical approximations that FLAME 
relies on are available for the full iJ-field formulation 
and involve Bessel / Hankel functions [20, 21, 22]. Two 
Bloch conditions - one for the H field and another one 
for the E field - need to be imposed on the cell bound- 
aries; the implementation details are described in [24]. In 
matrix-vector form, the FLAME eigenvalue problem is 
LE_ = [b^Bx + byBy)E_, where E_ is the Euclidean vector 
of the nodal values of the field on the grid, and b^, by are 
the Bloch factors b^ — exp{iKxL.x); by — exp{iKyLy). 
Matrix L has a sparse structure typical of finite-difference 
methods; matrix B is extremely sparse - its nonzero en- 
tries correspond only to two layers of grid nodes adjacent 
to the cell boundary. 

Numerical examples for real-K modes in regular di- 
electrics are given in [24]. To illustrate the compu- 
tation of complex modes in plasmonic crystals, let us 
start, for convenience of comparison, with the same ex- 
ample as in [4]. The permittivity of the rods is as- 
sumed to be described by the normalized Drude model 
e((D) = 1 — uj^^{lu — id), 

CieS OJ ^ LOfuJpj UJc -- ^^c/^^p^ vviicic ijjp 

frequency, lOc is the damping rate, and c is the speed of 
light in free space. The square lattice cell size is a = c/ojp 
and the cylinder radii are riod = 0.45a. The Bloch modes 
computed with the 9-point FLAME scheme coincide with 
the ones reported in [4] ; a fragment of the band diagram 
computed with FLAME for Wc = is shown in Fig. 1 for 
reference. The accuracy of FLAME is remarkably high: 
as evidenced by the typical convergence results for sev- 
eral modes in Table I, on the 40 x 40 grid the FLAME 



I ^ with the normalized frequen- 
ce /wp, where tUp is the plasma 
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FIG. 2: A backward-wave mode (TE3, circles, dotted line) 
in a plasmonic crystal with the lattice cell size as small as 
~ Ao/40. Parameters: ujp — 0.25, rrod ~ 0.45a. 
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Krcsl 


KC 


K^ 


20 X 20 


0.241048852 


2.573701526 


4.113530131 


30 X 30 


0.240708525 


2.572916824 


4.108636747 


40 X 40 


0.240652136 


2.572740104 


4.544368769 


50 X 50 


0.240644493 


2.572711175 


4.544040773 



TABLE I: Typical convergence of the numerical values for 
Bloch wavenumbers. Results for a) = 0.26. A'rcai is a real 
mode. Modes /s'1'2 have the real part K' = n/a. 



results already have about five correct digits - several 
orders of magnitude higher accuracy than in a typical 
PWE calculation [11, 22, 24]. Note that this example is 
not computationally favorable because the filling factor 
is high and the gaps between the rods are narrow. 

There is some similarity between FLAME and the 
semi- analytical construction of [28], where the field in 
the whole lattice cell is expanded into cylindrical har- 
monics and an eigenvalue problem for the coefficients of 
this expansion is obtained by imposing the Bloch bound- 
ary conditions at a set of collocation points on the cell 
boundary. In contrast, FLAME uses only local analyt- 
ical approximations that are valid over each individual 
grid stencil. Therefore FLAME can be extended to more 
general situations, e.g. with more than one rod in the 
cell, non-circular shapes, defects, photonic waveguides, 
etc. [3, 10, 22], where global analytical solutions are not 
available. 

Backvifard waves. As shown in [23], the cell size of 
periodic dielectric structures capable of supporting back- 
ward waves must lie above a fundamental threshold: 

[ryl < (47r2)-i|e|„,,,(l + |A|„,ax(/:r')|e|max) (4) 

where r/ = o)^^; lo ~ uja/c = 27ra/Ao; |A|niax('C7^) is 
the maximum eigenvalue of the inverse of the electro- 
static operator £<; = V • eV. This eigenvalue is bounded 
unless the operating frequency is close to the quasi- 
static plasmon resonance value. In particular, for non- 
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FIG. 3: Predominantly negative power density on the right 
edge of the very small lattice cell (a/Ao — 0.02387). 



plasmonic materials with < Cmin < e < Cmax through- 
out the lattice cell, the lower cell size bound specializes 
to (a/Ao)2 > [|eUax{l + |e|max/(47r2e^i„))]-i [23]. The 
plasmonic case is thus exceptional, because at the plas- 
monic resonance the electrostatic operator in (4) becomes 
singular and the constraint on the cell size disappears. 

The following example demonstrates that backward 
waves, and hence negative refraction, can be obtained 



in plasmonic crystals with very small lattice cells (here 
^ Ao/40). As before, consider a square lattice of cylindri- 
cal rods, with rmd = 0.45a. For this "proof-of-concept" 
example, losses are neglected and the Drude-like dielec- 
tric function is set as Crod = 1 — <^p/a)^, with ujp — ujpa/c. 
A fragment of the band diagram in the TX direction 
for ujp — 0.25 is shown in Fig. 2. The TE3 band (circles, 
solid line) corresponds to a backward wave: (a) the group 
velocity dcu/dKB for this band is negative; (b) the x com- 
ponent of the Poynting vector, plotted in Fig. 3, clearly 
indicates negative energy flow; (c) the first-Brillouin-zone 
component of the wave is dominant, as can be confirmed 
by Fourier analysis of the periodic field i?pER(a;, y) in (2); 
hence the wave has a well defined positive phase velocity. 



In summary, negative refraction can be obtained in a 
quasi-homogeneous medium with a very small cell size 
as a "plasmonic exception" circumventing the theoreti- 
cal cell size bounds of [23]. This may have important 
implications for imaging, as the inhomogeneity of the 
medium limits the resolution [15, 23]. For illustration, 
a backward-wave mode is demonstrated in a plasmonic 
crystal with the cell size of less than l/40th of the vac- 
uum wavelength. FLAME, a generalized finite-difference 
calculus based on very accurate local approximations of 
the solution, was applied to the computation of real and 
complex Bloch bands. 
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